#ifndef __SPARSEMATRIX__H__

#define __SPARSEMATRIX__H__

void sparse_add(int **J0, double **S0, int row, int col, double val);

void sparse_compress(int **J0, double **S0, int *I, int *J, double *S, int dim);

double sparse_entry(int i, int j, int *I, int *J, double *S, int dim);

void GS(int *I, int *J, double *S, double *rhs, double *x, int dim);

int CG(int *I, int *J, double *S, double *b, double *x, int dim);

#endif
/**
 * @file   sparse_matrix.c
 * @author Xian-Liang HU
 * @date   Tue Jan  1 10:32:03 2008
 * 
 * @brief  稀疏矩阵的存储以及迭代求解方法子程序
 * 都是自己写的哦，顺便熟练一下稀疏线性方程组的求解，
 * 存储的格式是通常用的行为主的存储方式
 */



//////////////////////////////// sparse matrix operation /////////////////////
/**
 * Sparse matrix add operation. The sparse matrix is represent by matrix J
 * and S.  Where J is the column index whose size is dim x MAXCOUPLE, and
 * the entries are double type and it's meaning is the entry of represented
 * Matrix. We assume that the vale of MAXCOUPLE is big enough and the first
 * column of J represent the number of nonzeros in current row.  so the
 * first column is initialized to be 1, it says that the pointer is at the
 * first position.
 */
void sparse_add(int **J0, double **S0, int row, int col, double val)
{
  int pos;
  // find if value col is in the row's row of matrix J
  for(pos=1; pos <= J0[row][0]; pos++){
    if(J0[row][pos] == col) break;
  }
  // whether J exist or not, i is the current position to insert
  // since S is initialized to be zero, add on is safe, also the rewrite of
  // J[row][i] is also safe!! This implementation is wonderful!
  J0[row][pos] = col;
  S0[row][pos] += val;
  J0[row][0] = (pos > J0[row][0])?pos:J0[row][0];
}

/**
 * treat the sparsity pattern from matrix J0 and S0 the row compress
 * storage format is gotten : I,J,S we assume the memory is reserved by
 * calling functions.  for ease of matrix iteration method, the diagnal
 * entries are put at the beginning of each row.
 */
void sparse_compress(int **J0, double **S0, int *I, int *J, double *S, int dim)
{
  int i,j,pos;
  // fill vector I first, whose first element must be 0
  I[0] = 0;
  for(i = 1; i <= dim; i++){
    I[i] = I[i-1] + J0[i-1][0];
  }
  // then run over J0 and S0 to fill J and S
  for(i = 0; i < dim; i++){
    for(j = 1,pos = 1; j <= J0[i][0]; j++){
      // remember that pos=0 is for the diagnal element
      if(J0[i][j] == i){
	J[I[i]] = i;
	S[I[i]] = S0[i][j];
      }
      else{
	J[I[i]+pos] = J0[i][j];
	S[I[i]+pos] = S0[i][j];
	pos++;
      }
    } 
  }
}

/**
 * get the specific entry value on position (i,j) 
 */
double sparse_entry(int i, int j, int *I, int *J, double *S, int dim)
{
  // for ease, we assume that i and j are not beyond dim
  // but for purpose of later use, we place dim here
  int pos;
  if (i==j){ // the diagnal element is easy to get!
    pos = I[i];
    return S[pos];
  }
  // otherwise the other element
  for(pos = I[i]+1; pos < I[i+1]; pos++)
    if(J[pos] == j) break;
  if(pos == I[i+1])
    return 0.0; // (i,j) are not int the pattern, so just return 0.0
  else
    return S[pos];  
}

/**
 * Solving linear system using various sparse matrix solver.
 * Currently, gauss-sediel works well for elliptic problem.
 */
void GS(int *I, int *J, double *S, double *rhs, double *x, int dim){
  double err = 1.0,res;
  int iter = 1,i,j,l;
  // double *x_old = (double *) calloc(dim,sizeof(double));
  while(err > 1e-12){
    // this is gauss sediel iteration:
    for(i = 0; i < dim; i++){
      res = rhs[i];
      for(l = I[i] + 1; l < I[i+1]; l++)
	res -= x[J[l]]*S[l];
      x[i] = res/S[I[i]];
    }
 
    //calculate the norm of the residual
    err = 0.0;
    for(i = 0; i < dim; i++){
      res = rhs[i];
      for(l = I[i]; l < I[i+1]; l++)
	res -= x[J[l]]*S[l];
      err += res*res;
    }
    err = sqrt(err);

    //if((iter++ % 100) == 0)
    //	printf("=== Iteration %d: the error is %le ===\n",iter-1, err);
  }
  for(i=0;i<dim;i++){
    printf("%le.\n",x[i]);}
  // printf("The G-S iterative solver finished at step %d with residual = %le.\n",iter-1, err);
}

/* 
 * CG solves the symmetric positive definite linear
 */
int CG(int *I, int *J, double *S, double *b, double *x, int dim)
{
  int max_iter,i,k,l;
  double *p,*r,*q,tol,alpha,beta;
  double nr,nr0,npq,val;

  max_iter = 3000;
  tol = 1e-12;
  
  p = (double*) calloc(dim,sizeof(double));
  q = (double*) calloc(dim,sizeof(double));
  r = (double*) calloc(dim,sizeof(double));
  // set r = b - A*x
  for(i = 0; i < dim; i++){
	r[i] = b[i];
	for(l = I[i]; l < I[i+1]; l++)
		r[i] -= S[l]*x[J[l]];
  }
  // calculate r0 and prepare p1
  nr0 = 0.0;
  for(i=0;i<dim;i++){
	nr0 += r[i]*r[i];
	p[i] = r[i];
  }
  if (nr0 <= tol){
	free(p);
	free(q);
	free(r);
	return 0; // convergence without iteration
  }
  
  // begin iteration
  for (k = 1; k <= max_iter; k++) {    
    // q = p*A*p
    npq = 0.0;
    for(i = 0; i < dim; i++){
	q[i] = 0.0;
	for(l = I[i]; l < I[i+1]; l++)
	    q[i] += S[l]*p[J[l]];
	npq += p[i]*q[i];
    }	
    alpha = nr0 / npq;
    // update x and r
    for(i=0;i<dim;i++){
	x[i] += alpha * p[i];
        r[i] -= alpha * q[i];
    }
    // calculate new residual
    nr = 0.0;
    for(i=0;i<dim;i++)
		nr += r[i]*r[i];
    printf("Iteration %d: residual = %le, alpha = %lf.\n",k,nr,alpha);
	
    // test convergence
    if (nr <= tol){
	printf("I am convergenced!\n");
	free(p);
	free(q);
	free(r);
        return k; // return the iteration count when convergence
    }
    // if not convergence, iteration again!
    beta = nr/nr0;
    for(i=0;i<dim;i++)
	p[i] = r[i] + beta * p[i]; // update p  
    nr0 = nr;
  }
  // otherwise, not convergence, return negative value
  free(p);
  free(q);
  free(r);
  return -max_iter;	
}
